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SUMMARY 

A  continuum  sensitivity  analysis  is  presented  for  large  inelastic  deformations  and  metal  forming  processes. 
The  formulation  is  based  on  the  differentiation  of  the  governing  field  equations  of  the  direct  problem  and 
development  of  weak  fonns  for  the  corresponding  field  sensitivity  equations.  Special  attention  is  given  to 
modelling  of  the  sensitivity  boundary  conditions  that  result  due  to  frictional  contact  between  the  die  and 
the  workpiece.  The  contact  problem  in  the  direct  deformation  analysis  is  modelled  using  an  augmented 
Lagrangian  formulation.  To  avoid  issues  of  non-differentiability  of  the  contact  conditions,  appropriate  regular¬ 
izing  assumptions  are  introduced  for  the  calculation  of  the  sensitivity  of  the  contact  tractions.  The  proposed 
analysis  is  used  for  the  calculation  of  sensitivity  fields  with  respect  to  various  process  parameters  including 
the  die  surface.  The  accuracy  and  effectiveness  of  the  proposed  method  are  demonstrated  with  a  number  of 
representative  example  problems.  In  the  die  design  applications,  a  Bezier  representation  of  the  die  curve  is 
introduced.  The  control  points  of  the  Bezier  curve  are  used  as  the  design  parameters.  Comparison  of  the 
computed  sensitivity  results  with  those  obtained  using  the  direct  analysis  for  two  nearby  dies  and  a  finite 
difference  approximation  indicate  a  very  high  accuracy  of  the  proposed  analysis.  The  method  is  applied  to 
the  design  of  extrusion  dies  that  minimize  the  standard  deviation  of  the  material  state  in  the  final  product 
or  minimize  the  required  extrusion  force  for  a  given  reduction  ratio.  An  open-forging  die  is  also  designed 
which  for  a  specified  stroke  and  initial  workpiece  produces  a  final  product  of  desired  shape.  Copyright 
©  2000  John  Wiley  &  Sons,  Ltd. 
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1.  INTRODUCTION 

Many  metal  alloys  are  currently  employed  in  the  manufacture  of  automotive,  aerospace  and  other 
hardware  components.  The  high  cost  of  manufacturing  critical  structural  components  can  be  sig¬ 
nificantly  reduced  with  the  development  of  mathematically  and  physically  sound  computational 
methodologies  for  process  design  and  control.  The  complicated  nature  of  polycrystalline  materials 
and  the  induced  changes  in  their  microstructure  during  processing  are  among  the  main  challenges 
that  one  must  consider  while  developing  means  for  the  design  and  control  of  bulk  forming  pro¬ 
cesses  that  result  in  products  of  desired  shape  and  material  state. 

The  desired  objectives  for  a  single  forming  operation,  e.g.  in  a  hot  extrusion  process,  may 
include  one  or  more  of  the  following  criteria:  uniform  deformation  in  the  final  product;  minimum 
required  work  or  extrusion  pressure;  desired  microstructure  in  the  final  product;  minimum  or  desired 
residual  stress  distribution;  minimum  deformation  and  wear  of  the  die;  desired  shape  of  the  final 
product;  and  minimum  porosity  in  the  final  product.  Any  of  these  objectives  can  theoretically  be 
achieved  by  appropriate  design  of  the  die  surface;  design  of  the  preform;  design  of  the  material 
state  (microstructure)  in  the  initial  billet;  and  appropriate  selection  of  the  process  parameters  (ram 
speed  and  pressure  history,  operating  temperature,  etc.).  However,  it  is  important  to  note  that  in  a 
single  forming  operation  there  is  only  a  limited  control  of  the  material  state  in  the  final  product 
that  one  can  achieve  using  a  single  stage  design  and  generally  a  multistage  process  design  is 
required. 

Most  deformation  process  design  is  currently  focussed  on  trial  and  error  techniques  based  on 
previous  experience  and  the  results  of  the  direct  analysis.  A  systematic  review  of  such  problems 
is  given  in  Reference  [1],  The  design  approaches  reported  in  this  reference  are  not  mathematically 
rigorous  and/or  realistic  from  a  material  representation  point  of  view.  However,  a  variety  of 
important  forming  design  problems  were  addressed. 

On  the  other  hand,  sensitivity  analysis  and  optimization  theory,  provide  a  fresh  look  at  these 
problems  and  can  lead  to  realistic  and  accurate  designs.  Traditional  forming  design  problems  that 
can  be  analysed  as  optimization  problems  are  the  optimum  design  of  dies,  preforms  and  process 
parameters.  To  mathematically  address  such  problems,  one  needs  to  calculate  the  sensitivity  of  the 
material  state  and  geometry  at  various  stages  of  deformation  with  respect  to  infinitesimal  changes 
in  each  of  the  design  variables. 

Sensitivity  computations  can  be  performed  using  finite  difference  approximations  and  the  results 
of  the  direct  analysis  for  two  nearby  design  variables.  In  addition  to  significant  computer  resources 
required  for  solving  the  direct  problem  multiple  times,  difficulties  arise  in  such  calculations  from  the 
fact  that  many  direct  analysis  tools  are  insensitive  to  infinitesimal  changes  in  the  design  variables 
(e.g.  the  die  surface)  and  cannot  provide  accurate  sensitivity  fields.  This  is  particularly  true  when 
the  computed  sensitivity  fields  are  of  the  same  order  of  magnitude  as  the  numerical  error  in  the 
solution  of  the  direct  analysis. 

In  the  direct  differentiation  method,  a  set  of  field  equations  are  developed  by  considering  the 
variation  of  the  continuum  or  discretized  field  equations  of  the  direct  problem  with  respect  to 
small  changes  in  design  parameters  [2-6].  The  sensitivity  field  equations  are  linear  and  can  be 
efficiently  solved.  Direct  automatic  differentiation  can  also  be  applied  to  evaluate  sensitivities  [7], 

In  most  of  the  current  developments  in  the  field  of  sensitivity  analysis  for  large  deformation 
problems,  the  flow  formulation  is  utilized  and  emphasis  is  given  to  steady-state  forming  applications 
[6].  Isothermal  as  well  as  non-isothermal  deformations  have  been  considered  and  the  sensitivity 
of  various  forming  processes  with  respect  to  shape  and  material  parameters  has  been  studied 
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[8—10] .  Computation  of  sensitivity  coefficients  for  the  purpose  of  parameter  identification  for  finite 
elasto-plastic  deformations  is  reported  in  Reference  [11], 

Although  the  last  few  years  have  seen  an  increasing  interest  in  the  direct  differentiation  method 
for  calculating  sensitivities,  a  number  of  challenges  still  remain  in  the  accurate  and  effective 
computation  of  sensitivity  fields  for  practical  engineering  problems  that  involve  frictional  contact 
[12],  A  sensitivity  analysis  for  frictional  metal  forming  processes  in  steady  state  using  a  flow 
formulation  is  presented  in  Reference  [8]  where  the  effect  of  variation  in  the  Coulomb  friction 
coefficient  on  the  deformation  response  is  studied.  Special  contact  elements  were  introduced  in 
the  workpiece  boundary  to  handle  contact/frictional  effects.  More  recently,  sensitivity  analysis  of  a 
non-steady  state  open  die  forging  process  has  been  carried  out  in  Reference  [13]  to  study  the  effect 
of  variation  in  the  height  of  the  initial  cylindrical  preform  in  a  single  stage  operation.  A  shape 
optimization  scheme  for  non-steady  state  metal  forming  processes  with  applications  to  preform 
design  in  forging  is  given  in  References  [14,  15].  Applications  of  optimization  techniques  to  net- 
shape  manufacturing  involving  forging  processes  are  given  in  Reference  [16].  Other  applications 
to  extrusion  and  forging  die  design  can  be  found  in  References  [17,  18]. 

The  main  objective  of  this  paper  is  to  provide  sensitivity  field  calculations  that  can  be  used 
for  accurate  forming  process  design.  A  Lagrangian  sensitivity  analysis  is  developed  to  include  the 
effects  of  contact  and  friction  for  both  steady  state  and  non-steady-state  metal  forming  processes. 
Emphasis  in  this  paper  is  given  to  non-shape  forming  design  problems  such  as  the  design  of 
various  process  parameters.  In  such  problems,  for  example,  one  calculates  the  die  surface  (for  a 
given  reduction  or  die  stroke)  that  results  in  desired  material  properties  in  the  final  product  or  a 
product  of  desired  shape  (e.g.  in  an  open-die  forging  process). 

It  is  common  to  discretize  the  die  surface  using  basis  functions  (e.g.  splines  or  polynomial 
functions)  assuming  that  certain  constraints  are  satisfied  (e.g.  fixed  reduction)  [19-22],  With  such 
discretizations,  a  die  design  problem  can  be  posed  as  an  optimization  problem  with  respect  to 
a  finite  number  of  algebraic  design  parameters.  The  design  parameters  are  here  denoted  by  the 
vector  P p  —  {Pi},  i  =  l,2,...  ,NP.  Various  problem-dependent  objective  functions  are  usually  written 
implicitly  in  terms  of  the  parameters  P/;.  Such  objective  functions  may  represent  for  example,  the 
deviation  of  the  resulting  material  state  of  the  final  product  from  a  desired  state  for  the  process 
defined  by  design  parameters  \\p.  Similarly,  the  objective  function  can  be  selected  to  minimize  the 
required  work  expenditure  for  the  design  parameters  P/;. 

There  are  many  ways  to  solve  the  above  class  of  optimization  problems  but  most  commonly 
a  sequential  search  method  is  employed  starting  from  a  guess  solution  that  is  iteratively  updated 
along  some  specific  direction  and  a  given  step  size  [23],  To  determine  at  each  optimization  iteration 
the  specific  search  direction  and  step,  one  must  evaluate  the  gradient  of  various  fields  (material  or 
geometric)  involved  in  the  definition  of  the  objective  function  with  respect  to  design  parameters. 
The  interest  here  is  in  using  a  sensitivity  analysis  for  evaluating  the  gradient  of  the  objective 
function.  With  such  an  analysis,  sensitivities  of  various  Lagrangian  fields  can  easily  be  calculated 
sequentially  in  time. 

The  work  reported  here  is  an  extension  of  an  earlier  work  by  Zabaras  and  co-workers  [24],  In  this 
work,  a  continuum  sensitivity  analysis  for  hyperelastic  viscoplastic  deformations  was  introduced 
and  some  die  design  examples  were  analysed.  The  sensitivity  fields  were  defined  as  directional 
derivatives  of  the  corresponding  Lagrangian  fields  of  the  direct  analysis.  For  die  design,  for  ex¬ 
ample,  these  directional  fields  were  calculated  at  a  reference  die  and  for  a  given  die  perturbation. 
Appropriate  kinematic  and  constitutive  sensitivity  problems  were  defined  and  the  sensitivities  of 
the  contact  tractions  were  calculated  by  assuming  Coulomb  friction  in  the  contact  boundary  and 
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by  applying  direct  differentiation  of  the  impenetrability  constraint  and  contact  traction  components 
with  respect  to  the  die  surface.  The  sensitivity  of  the  contact  pressure  was  eliminated  from  these 
equations  with  the  sensitivity  of  the  normal  to  the  die  unit  vector  expressed  in  terms  of  the  refer¬ 
ence  die  surface  and  the  die  surface  variations.  Even  though  the  technique  introduced  in  Reference 
[24]  was  found  to  give  very  good  results  for  simple  die  design  problems,  accuracy  concerns  were 
raised  in  many  other  die  and  preform  design  problems  [25],  These  difficulties  were  attributed  to 
the  explicit  manner  with  which  contact  was  treated  in  both  the  direct  and  sensitivity  analyses. 
Such  issues  will  be  addressed  in  this  paper  by  employing  a  fully  implicit  contact  scheme  and  its 
regularized  derivative. 

The  plan  of  this  paper  is  the  following.  At  first,  a  brief  introduction  of  the  direct  defor¬ 
mation  problem  is  given  to  set  the  notation  for  the  sensitivity  analysis.  A  fully  implicit  aug¬ 
mented  Lagrangian  algorithm  for  handling  the  contact/friction  sub-problem  of  the  direct  analysis 
is  then  presented.  The  constitutive  and  kinematic  sensitivity  problems  are  briefly  reviewed  next 
together  with  the  techniques  for  their  solution.  The  contact  sensitivity  algorithm  is  derived  by  the 
direct  differentiation  of  continuum  contact  and  frictional  constraints  and  the  implementation  of 
the  resulting  sensitivity  constraints  using  a  penalty  formulation.  To  allow  for  the  differentiability  of 
the  continuum  contact/frictional  constraints,  appropriate  regularizing  assumptions  are  introduced. 
The  developed  contact  sensitivity  problem  is  linear  and  augmentation  steps  are  avoided  by  ap¬ 
propriate  selection  of  penalty  parameters  that  enforce  contact/frictional  sensitivity  constraints.  The 
paper  concludes  with  an  evaluation  of  the  accuracy  of  the  proposed  methodology  and  a  presentation 
of  the  solution  of  a  number  of  die  design  problems  in  extrusion  and  open-die  forging  processes. 


2.  REVIEW  OF  A  DIRECT  LARGE  INELASTIC  DEFORMATION  PROBLEM 
2.1.  Kinematic  and  constitutive  equations 

The  basic  governing  equations  for  the  Lagrangian  analysis  of  the  large  deformation  of  isotropic 
viscoplastic  solids  are  presented  here  to  set  the  notation  for  the  proposed  sensitivity  analysis. 

Let  us  denote  with  B0  the  configuration  of  the  body  at  time  t  =  0  and  with  B,  the  current  body 
configuration  at  time  t.  Consider  the  location  x,  in  the  current  configuration  B;,  that  at  time  t  is 
occupied  by  the  material  point  p.  Let  X  be  the  location  of  this  particle  in  the  configuration  B0  at 
time  t  =  0.  A  smooth  mapping  cj)  exists  such  that 

x  =  4>(X,o  (i) 

The  deformation  gradient  F  is  then  defined  by 

F  =  V<KX,f)  (2) 

Assuming  isothermal  conditions,  the  deformation  gradient  F  is  decomposed  as  follows: 

F  =  FeFp,  det  Fe  >  0  (3) 

where  Fe  is  the  elastic  deformation  gradient  and  Fp  the  plastic  deformation  gradient  with  det  Fp  =  1 . 
From  the  polar  decomposition  of  Fe,  one  can  write, 

Fe  =  ReUe  (4) 


Copyright  ©  2000  John  Wiley  &  Sons,  Ltd. 


hit.  J.  Numer.  Meth.  Engng  2000;  48:679-720 


METAL  FORMING  PROCESSES  WITH  APPLICATIONS  TO  DIE  DESIGN  PROBLEMS 


683 


The  equilibrium  equations  can  be  expressed  as, 


V-P  +  f  =  0  VXeBo  (5) 

where  P(X,  l)  is  the  first  Piola-Kirchhoff  stress  and  f(X,  l)  is  the  body  force  density  in  the  reference 
configuration  Bo.  The  definition  of  P  and  f  is  as  follows: 

P(X,  t)  =  det  F  T  F_t 

(6) 

f(X,t)  =  detFb 


Here,  T  is  the  Cauchy  stress  and  b  the  known  body  force  density  in  the  current  configuration. 

The  weak  form  of  the  above  equilibrium  equations  is  stated  as  follows:  Calculate  <)>(X,  t)  such 
that 


*■£<**= 


PB„ 


t  •  udH,i+]  + 


b  •  u  d  V„ 


H+ 1 


(7) 


for  every  admissible  test  function  u.  The  hyperelastic  constitutive  equations  are  taken  as  [26], 

T  =  ^e(Ee)  (8) 

where  the  elastic  strain  Ee  is  defined  by 

Ee  =  ln(Ue)  (9) 

and  the  corresponding  conjugate  stress  T  is  given  as, 

T  =  (det  Ue )  (Re  )T  TRe  (10) 

The  isotropic  elastic  moduli  are  defined  by 


=  +  (jt- 


(ii) 


with  :'/f  the  bulk  modulus,  Ct  the  shear  modulus  and  I,  ■/,  the  unit  second-  and  fourth-order 
tensors,  respectively. 

A  flow  rule  is  given  in  the  form  of  the  evolution  of  Fp  with  zero  spin  of  the  intermediate 
configuration  [26], 


where 


Lp=  Dp  =  Ep(Fp)_1  = 


^  ep  Np(  T',  5) 


NP(T',  a) 


the  equivalent  stress  a  is  given  by 


with 


(12) 


(13) 


(14) 


(15) 
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Finally,  the  equivalent  plastic  strain  rate  ;TP  is  specified  as 

£p  =  /(<7,5)  (16) 

and  the  evolution  of  the  scalar  state  variable  s  (internal  resistance  to  plastic  deformation)  is  given 
by 


s  =  g{5,s)  =  h{5,s)f{5,s)  (17) 

The  functions  /(<r,s)  and  h{a,s)  are  experimentally  determined  for  a  particular  material  [27], 
With  the  above  constitutive  equations,  one  can  solve  the  weak  form  of  the  equilibrium  equations 
(Equation  (7))  to  calculate  the  deformation  and  evolution  of  state  induced  by  given  applied  loads 
or  kinematic  boundary  conditions  [28].  Contact  constraints  play  an  important  role  for  metal  forming 
processes  [29,  30],  A  brief  presentation  of  the  contact  sub-problem  is  given  next  in  order  to  set 
the  notation  for  the  development  of  the  contact  sensitivity  algorithm. 


2.2.  The  contact  sub-problem 

A  schematic  of  the  contact  sub-problem  is  shown  in  Figure  1.  The  presentation  is  limited  to  plane- 
strain  and  axisymmetric  problems  and  dies  are  considered  to  be  rigid.  The  die  Q>  is  parametrized 
using  a  parameter  £  and  the  die  functions  y(£)  =  (yi(0, y3(0)^ 0<  <[;  <  1.  A  fixed  right-handed  basis 
(ei,e2,e3)  is  defined,  with  e2  pointing  into  the  plane  of  the  paper.  A  convected  basis  (r,e2,v)  is 
defined  at  each  point  (i.e.  for  each  particular  value  of  £).  The  tangent  vector  X\,  and  the  unit 
tangent  vector  r  are  given  as 


Sy i  ,  dy3 

t,=y,i  =  we|  +  — 


e3, 


(18) 


The  die  separates  the  space  into  admissible  and  inadmissible  regions  and  the  die  is  parametrized 
such  that  the  normal  vector  v  is  pointing  into  the  admissible  region  [29],  With  this  convention, 
the  unit  normal  vector  v  points  towards  the  body  (i.e.  v  =  — n,  where  n  is  the  outward  unit  normal 
vector  to  the  body).  Figure  1  shows  a  schematic  of  the  contact  sub-problem  and  introduces  the 
notation  for  the  various  contact  parameters.  The  gap  function  g  of  any  point  x  in  space  is  defined 
as  the  shortest  distance  of  that  point  from  the  die.  Thus,  we  write 


y  —  x  =  gf(x)v(y)  (19) 

where  y  £  Q)  is  the  value  of  y  that  minimizes  the  norm,  ||x  — y||.  A  unique  value  of  the  parameter  f 
is  associated  with  each  y. 

Following  the  work  in  Reference  [29],  the  contact  traction  vector  X  per  unit  reference  area 
r  c  3B„  is  introduced  and  its  components  zN  and  zT  are  defined  as  follows: 


k  =  7,NV  —  (20) 

The  impenetrability  constraints  can  now  be  written  as  follows:  For  all  x„  £  I  ,  with  x„+i  = 

x„  +  u(x„), 


2n^0,  gf(x„+i)sC0,  2Nq(x„+i )  =  0 


(21) 
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Admbrihlc  K  (  jl  *  fl  J 

(■' 


Figure  1.  A  schematic  of  the  contact  sub-problem  showing  the  body  configurations  at  times  t„  and  t„+ 1,  the 
definition  of  the  gap  function  g(x)  and  the  admissible  and  inadmissible  regions.  I’  refers  to  the  part  of  the 

boundary  8B„  that  could  potentially  come  in  contact. 


The  Coulomb  friction  law  can  be  written  as: 

Lj  =  —  L  +  XnV 
T  :=  ||^t||  —  /(^-n  C  0 

VT  =  (22) 

1 1 A/T II 

i  >  0 

ZT  =  0 

An  augmented  Lagrangian  algorithm  that  enforces  the  above  contact  and  friction  constraints,  is 
presented  in  the  next  section. 

2.3.  Basic  elements  of  a  Lagrangian  analysis 

In  the  updated  Lagrangian  FEM  formulation,  a  sequence  of  incremental  problems  is  defined  from 
time  t„  to  4+i,  n  =  0, 1,2, ...  .  Within  each  time  increment,  the  deformation  problem  is  further  di¬ 
vided  into  three  incremental  sub-problems:  (1)  the  constitutive  problem,  (2)  the  kinematic  problem 
and  (3)  the  contact  problem. 

In  the  constitutive  incremental  problem,  one  calculates  (T„+i, ,4+1,1^,  , )  given  the  incremental 
deformation  gradient  Fr  from  t„  to  4+1  (see  Figure  1)  and  (T(I,,s'„,  F„p).  Constitutive  integration 
techniques  for  such  problems  can  be  found  in  Reference  [31].  In  the  present  analysis  the  radial 
return  mapping  of  Reference  [26]  is  implemented. 
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In  the  kinematic  problem,  one  calculates  the  incremental  displacements  from  tn  to  tn+\  given 
the  triad  (T,s,Fp)  at  tn+\.  A  Newton-Raphson  scheme  is  developed  by  linearizing  the  weak  form 
of  the  equilibrium  equation  [28]. 

In  the  contact  sub-problem,  given  the  configuration  B„+i,  the  die  location  and  shape  at  tn+ as 
well  as  estimates  of  the  contact  tractions  from  a  previous  time  step  or  a  Newton  iteration,  one 
computes  (updates)  regions  of  contact  as  well  as  the  contact  traction  components  zN  and  /,+. 

The  time  integration  of  the  contact  sub-problem  based  on  the  augmented  Lagrangian  analysis 
of  Laursen  and  Simo  [29,  30],  is  briefly  reviewed  since  it  is  essential  for  the  presentation  of  the 
proposed  sensitivity  methodology.  The  time  integration  of  the  frictional  constraint  is  achieved  by 
the  introduction  of  a  trial  state  and  a  subsequent  return  map.  In  the  algorithm  below,  k  refers  to 
the  augmentation  index  and  j  refers  to  the  Newton-Raphson  iteration  index.  The  terms  G  and  Gc 
in  the  kinematic  problem  refer  to  the  principle  of  virtual  work  and  contributions  from  the  internal 
work/non-contact  related  boundary  terms  and  contact  terms,  respectively  [28], 

1.  Initialization: 

4°)  =  (^n„  +£n0(x^°+)1)> 

A2(t0)  =  0 
k  =  0 

2.  Solve  for  the  incremental  displacement  =  —  x„  which  is  the  converged  estimate  of 

u(*)_x(*)  x 

G(uf},  u)  =  G{uf\ u)  +  Gc(uf\ u)  =  0 

Note:  The  above  equation  is  non-linear  and  is  solved  iteratively  (iteration  j)  [28].  A  Newton- 
Raphson  algorithm  is  used  to  solve  this  system  and  the  computation  of  the  linearized  stiffness 
matrix  and  residual  is  dependent  on  the  solution  of  the  kinematic,  constitutive  and  contact  sub¬ 
problems.  The  solution  of  this  non-linear  equation  predicts  estimates  of  the  body  configuration 
Bf„+])  ■  at  4+i  and  finally  the  converged  solution  B*+1. 

The  contact  traction  ~kk-_ ,  used  in  the  solution  of  the  linearized  principle  of  virtual  work 
(based  on  the  estimate  B*b+1)  -_x)  is  calculated  as  follows: 

Normal  traction: 

An  =  (4T  +  £n 

Tangential  traction: 

Compute  trial  state: 

4'al  =  K  +  ST (l((„%- ,  -  In)  + 

Ttrial=||k![ial||  -/dN 

Radial  return  update: 

IF  (Ttnal^0) 

THEN  =  2+'al  (stick) 
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ELSE  (return  map) 

^  trial 

*t  =  A^n  ppjj  (slip) 

Complete  contact  traction  description: 
4-i  =  2Nv(y)  -  ATTi(y) 


3.  Check  for  contact  constraints  satisfaction: 

IF  gCx^XTOLg  AND 

ul+i  -  ln\\  ^TOLf  Vx„  G  r  such  that  ||kT||  </t(4 *  +  Sn^x^)} 

THEN 

CONVERGED  (GOTO  4) 

ELSE  AUGMENT  (Vx„  ef) 


4i+1)=(4i)+^(xii+)1)) 


y' trial 


4ial|l  -  //;.<*• 0 


IF  (Ttrial^0) 
THEN 


A4+1)  =  A4i}  +  £tA|W  (stick) 


ELSE  (return  map) 

}  trial 

a4+1)  =  ^4+I)  (slip) 

ENDIF 
k  =  k  +  1 
GOTO  2 
ENDIF 

4.  Post-process  operation: 

Use  the  converged  solution  of  the  configuration  B„+i  to  update  the  triad  (T,  s,  Fp)  as  well  as 
the  contact  traction  components  (2n,2t)  at  time  tn+ 


The  three  sub-problems  (constitutive,  kinematic  and  contact)  are  solved  in  an  iterative  manner. 
More  details  on  the  implementation  of  these  problems  can  be  found  in  References  [28,  32]  where 
an  object-oriented  environment  for  finite  deformation  problems  is  discussed. 


3.  THE  SENSITIVITY  PROBLEM 
3.1.  Definition  of  sensitivity  fields 

A  definition  is  given  for  the  sensitivity  of  any  field  related  to  the  direct  problem  (e.g.  of  the 
deformation  gradient,  stress,  state  variables,  etc.)  with  respect  to  the  die  shape.  The  parameter 
|5 p  introduced  below  is  either  the  die  surface  itself  (i.e.  P/(  is  a  function)  or  a  set  of  algebraic 
parameters  that  define  the  die  surface  in  a  finite  dimensional  space  as  discussed  in  Section  1. 
Before  we  proceed,  it  is  mentioned  that  §p  can  also  represent  any  other  non-shape  related  process 
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Updated  reference 


Figure  2.  Schematic  diagram  of  the  deformation  history  of  a  workpiece  corresponding  to  two  dies  defined 
by  the  design  parameters  |5/;  and  (|5,,  +  A|5/; ).  Note  that  in  the  deformation  process  with  a  die  [fr  +  Ap/;,  the 

relative  deformation  gradient  from  tn  to  tn+ 1  is  easily  shown  to  be  Fr+  Fr,  where  Fr  =(F«+i  —  Fr  F„)F“1. 


parameter  such  as  the  ram  or  die  speed,  operating  temperature,  lubrication  conditions  as  well  as  the 
initial  material  state.  For  preform  design  and  related  shape  optimization  problems,  see  References 
[25,  33],  For  clarity  of  the  presentation,  unless  otherwise  stated,  we  only  refer  to  sensitivities  with 
respect  to  the  die  surface. 

Let  <1>(X,  t;  Pp)  be  a  Lagrangian  field  defined  in  the  direct  analysis  for  a  given  choice  of  design 
P p.  The  notation  0(X,  f;P,,)  implies  that  can  be  uniquely  determined  by  solving  the  direct 
deformation  problem  for  a  particular  set  of  boundary  conditions  resulting  from  a  die  defined  by 
P p.  Figure  2  shows  the  configurations  of  the  body  during  a  deformation  process  for  two  different 

but  nearby  sets  of  design  parameters  (i.e.  two  nearby  dies). 

* 

We  define  the  sensitivity  <D  (X,  1;  P;),  Ap/; )  of  <D  as  the  directional  derivative  (Gateaux  differen¬ 
tial)  of  <1>  with  respect  to  Pp  in  the  direction  of  App,  i.e. 

<J>(X,f;  p^App)  =  A  [O(X,f;Pp  +  2Ap;,)]|/l  =  0 

=  0(X,f;  Pp  +  APp)  —  <1>(X,  t;  §p)  +  0(||  Ap^||2)  (23) 

Based  on  the  above  definition,  one  can  easily  define  temporal  and  spatial  derivatives  of  sensitivity 
fields  [24].  The  sensitivity  of  Eulerian  fields  can  be  calculated  by  first  transforming  them  to 
Lagrangian  fields  (via  the  mapping  given  by  Equation  (1))  and  then  applying  the  definition  given 
in  Equation  (23). 

3.2.  The  overall  deformation  sensitivity  problem 

Given  the  die  design  parameters  P;)  and  with  the  remaining  boundary  and  process  conditions  fixed, 
one  can  solve  the  direct  deformation  problem  and  calculate  the  material  state,  deformation  and 
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stresses  at  each  time  step.  On  the  other  hand  in  the  sensitivity  problem,  one  is  interested  in 
computing  the  sensitivity  history  of  the  material  state,  deformation  and  stresses  around  a  die  P;, 
and  with  a  die  perturbation  App.  In  simple  terms,  the  sensitivity  of  a  given  direct  field  tells  us  how 
much  this  field  will  change  from  its  value  obtained  from  the  direct  deformation  problem  solution 
obtained  with  a  die  §p  if  the  die  is  perturbed  by  an  amount  AP;). 

As  in  the  direct  analysis,  the  incremental  deformation  sensitivity  problem  (from  time  t„  to 
time  tn+ 1)  is  divided  into  three  coupled  sub-problems,  the  constitutive,  kinematic  and  contact 
sensitivity  sub-problems.  In  the  incremental  constitutive  sensitivity  problem,  given  the  sensitivity 
fields  at  time  tn,  [s„,F„],  as  well  as  the  direct  fields  at  time  tn+\,  [T„+1,s„+i  ,F„+i,F„p+1],  one 

*  *  *  * 

calculates  the  linear  relations  between  each  of  the  sensitivity  fields  [Th+i,  s«+i,Fjj+i]  and  F„+i-  In 
the  incremental  kinematic  sensitivity  problem,  given  the  direct  problem  solution  at  tn+\ ,  i.e.  given 

*  *  *  =t= 

[T„+i  ? F/j-|-i, Fn_j_ j ],  and  with  the  given  linear  relations  between  j  'IA  ■  i  -  V//  ■  i  -  If.'  ■  i  ]  and  F«+i, 
* 

one  calculates  F«+i-  In  the  incremental  contact  sensitivity  problem,  given  the  contact  traction 
components  In  and  At  at  time  tn+ \  as  well  as  the  identity  of  regions  of  sliding  and  sticking,  one 
computes  the  linear  relation  between  the  sensitivities  An  and  At  at  tn+\  and  the  sensitivity  x  of 
the  contact  boundary  in  the  current  configuration  B„+i. 

Remark  1.  For  simplification,  it  has  been  selected  in  the  direct  analysis  to  track  and  report  the 

calculation  of  the  triad  V„+i  =  [T„+i  ?  $  n+ 1 ?  Fp  J  and  F„+i,  whereas  in  the  sensitivity  analysis  the 
*  *  * 

duo  W„+|=  [s„+i,FJ+i]  and  F„+i  are  computed  and  reported.  The  selection  in  the  sensitivity 
*  * 
analysis  of  F®+i  to  define  the  sensitivity  of  the  intermediate  configuration  instead  of  FJ]+ 1  is  for 

*  *  - 

convenience  of  the  presentation  since  Fjj+i  and  F„+i  uniquely  define  F{+  for  given  Fp+1  and  F„+i. 


The  three  incremental  sensitivity  problems  (constitutive,  kinematic  and  contact)  are  linearly 
coupled  and  together  provide  a  linear  problem  for  the  calculation  of  sensitivities  of  both  the 
deformation  and  the  material  state  at  tn+\.  It  should  be  noted  that  the  corresponding  sub-problems 
in  the  direct  problem  are  non-linearly  coupled  and  are  solved  in  an  iterative  fashion  [31]. 


3.2.1.  The  constitutive  sensitivity  problem.  The  incremental  constitutive  sensitivity  problem  is 
defined  in  Reference  [24]  and  is  briefly  described  here  for  clarity  and  completeness  of  the  pre¬ 
sentation.  The  present  mathematical  approach  to  this  problem  is  similar  to  that  of  Reference  [24], 
however,  the  implementation  is  quite  different.  The  solution  of  the  constitutive  sensitivity  problem 
follows  the  following  steps: 

*  *  * 

Calculation  of  the  linear  relations  between  a„+\  and  FJ}+i 

* 

By  following  the  definitions  of  T'  (Equation  (15))  and  a  (Equation  (14)),  the  sensitivities  T[1+1 
* 

and  <f„+i  are  obtained  as  follows: 

Tn+i  =  T/j+i  —  -tr(T,i+i)I  (24) 


and 


*  _  3  T»+i  •  T«+i 

&n+ 1  0 

^  0/2+1 


(25) 
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* 

By  differentiating  the  hyperelastic  law  in  Equation  (8),  T„+i  is  expressed  as: 

t„+i  =  ^e[fc+i]  (26) 

* 

By  employing  the  first  Pade  approximation  of  the  logarithm  of  Ue  and  taking  the  sensitivity,  E®1+i 
can  be  written  as 

h+i  =  4(U„e+1  +  ir'U^CU^,  +  I)- 1  (27) 

where 

U„e+1  =  sym(U„e;ilsym(Ff+lF„e+, ))  (28) 

Equations  (24)-(28)  and  the  linear  transformation  equations  presented  in  Appendix  A  of  Reference 

*  * 

[24]  can  easily  be  combined  to  provide  the  linear  relations  between  T'n+\  and  F,®+i  as  well  as 
*  * 

<f„+i  and  FJj+1. 


-  _  1 

Calculation  of  the  linear  relation  between  Fj]+i[Fp]n+1  and  F^+[ 

According  to  the  definition  of  DP  (Equation  (12)),  one  can  write  the  following: 


DP  =  -^-(FP)[FP]_1  -  DPFp[Fpr1 
at 

=  -^-(FP[FP]_1)  +  FP[FP]_1  DP 
at 


DPFP[FP]_1 


(29) 


We  employ  an  Euler-backward  integration  scheme  over  (tn,tn+ 1)  to  integrate  Equation  (29)  and 
express  the  above  equation  as  follows: 


^(FPtFPJ-^  +  FP^]-1^^  -  DP+1FP[FP]_1  =DP+1 


Pre-multiplication  of  the  above  equation  by  exp(— tDp+1),  post-multiplication  by  exp(tDp+1)  and 
integration  of  the  resulting  equation,  leads  to  the  following: 


^+i[Fp]„“+i  =  AFPfStF^-^AFP)-1  +  At6P+1 

where  from  the  constitutive  sub-problem  of  the  direct  problem  [26], 

AFp=Fp+1[Fp]-1=exP(AtDp+1) 

with  At  =  tn+\  —  tn.  From  Equation  (12),  one  can  derive  the  following: 


f  * 
i«+l  rp/  , 
- -  1/1  + 1  ~T 


3 

2 


I  &n+l(f d)n+l  fn+l  ^  (fs)n+\  *  \  -fv 

[  ^2  an+l  +  -  5h+1  |  ^  //+1 

\  an+ 1  <7»+1  / 


(30) 

(31) 


(32) 
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where  from  Equation  (17),  s  is  obtained  from 

d  *  ^ 

—(s)  =  gss  +  gga  (33) 

at 

The  subscripts  to  the  functions  /(<t,s)  and  g(<J,s )  are  used  to  denote  partial  derivatives  with  respect 
to  their  arguments.  An  Euler-backward  integration  of  Equation  (33)  results  in  the  following: 

*  _  *  1  ,  (gf<f)n+l At  * 

5,1+1  “  5,1  1  -  A t(gs)n+l  +  1  -  At(p,)„+i  17,1+1 

Substitution  of  Equations  (32)  and  (34)  in  Equation  (30)  results  in 
Fp+i[Fp]„+11  =  C„+i  +  a„+iT,,,+[  +  bn+\dn+\i'n+{ 


(34) 

(35) 


where 


c„+1 = +  2,„+l(i(f'X..u,)  *'  A,t»« 

_  3/„+iA/ 

1  ^  ~ 

^  ®n-\- 1 

i  _  3  I  ®n+\(f  d^)n+\  fn+ 1  (/i)«+l  ($(t)h+1  ^  I  . 

"+1  2  |  <t„2+1  <f„+i(l  -  At(^)„+i)  J 


(36) 

(37) 

(38) 


The  sensitivity  fields  T'n+\  and  an+\  were  written  earlier  as  linear  functions  of  F,®+i  (see  Equa- 

* 

tions  (24)-(28)).  As  a  result,  Equation  (35)  provides  the  linear  relation  between  F}(+ 1 [Fp]~.'|  and 
* 

F„e+,. 


*  * 

Calculation  of  the  linear  relation  between  F,®+1  and  F«+i 

Starting  from  the  multiplicative  decomposition  of  the  deformation  gradient,  one  can  write 

F/,+i  =F„e+lFp+,  +F:+,Kp+,  (39) 

Hence, 


[FT+,(F«+1  [F];+,)F:+1  =  [Fe]^,F„e+l  +  |p+1[Fp] 


-1 
72+ 1 


(40) 


As  has  already  been  shown,  T'+i  and  an+\  linearly  depend  on  F^+) 

*  * 

Equation  (35)  in  Equation  (40)  results  in  a  linear  relation  between  F,®+i  and  F«+i 


Therefore  substitution  of 


[F  ]„+i(Fn+i  [F],1+1)Fn+|  C/;  1 1  [F  T  an^iTn^i  T  bfl-^ian^] Tw_j_^ 


(41) 
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Using  the  simplifications  of  Appendix  A  and  the  linear  transformation  equations  discussed  in  the 
Appendix  of  Reference  [24],  one  can  write  the  following: 

F„e+1  =^(V„+1)[F„+i]  +  A(V„+1,W„)  (42) 

where  A  is  a  second-order  tensor  function  and  si  a  fourth-order  tensor  function.  A  and  si  depend 
upon  the  direct  fields  \n+i  =  (T„+iAn+i,F^+1)  calculated  at  time  tn+\  and  the  sensitivity  fields 
*  *  * 

W„=  F,j)  calculated  at  time  tn. 


*  * 

Calculation  of  the  linear  relation  between  T«+i  and  FJj+i  Linearizations  similar  to  the  ones  per¬ 
formed  in  the  direct  analysis  as  part  of  the  calculation  of  the  tangent  stiffness,  are  used  here  to 

* 

calculate  the  sensitivity  of  the  Cauchy  stress.  Following  References  [24,  26],  one  can  express  T 
as 


* 

T«+i 


(det  Ue„+1)-1  R:+l  t„+,Re„T+1  -  tr(!e„+1)T,1+1  +  2sym(R?i+1Re„T+1T,1+1) 


(43) 


where 


R«+i[Re]I+i  =F^+1[Fe]„-+\  -  R^U^U^tM-H 


je-i-l 


el  —  1 


>enT 


(44) 


Expressions  for  E),+  |  and  U|+i  in  terms  of  F,j+1  are  given  in  Equations  (27)  and  (28),  respectively. 


The  calculation  of  T„+i  in  terms  of  F,j+i  was  reported  above  (Equations  (26)-(28)). 


3.2.2.  The  kinematic  sensitivity  problem.  To  calculate  the  sensitivity  field  F,  one  needs  to  consider 
the  sensitivity  of  the  equilibrium  Equation  (5)  in  conjunction  with  a  set  of  boundary  conditions 

for  this  kinematic  problem.  In  the  following,  a  principle  of  virtual  work  like  equation  is  developed 
* 

to  obtain  F-  The  non-contact  related  sensitivity  boundary  conditions  are  treated  in  this  section, 
whereas  the  contact  sensitivity  problem  is  addressed  in  the  next  section.  To  simplify  the  notation, 
we  omit  the  subscript  (n  +  1)  for  tensor  variables  that  refer  to  time  tn+\ . 

The  directional  derivative  of  the  equilibrium  equations  (Equation  (5))  results  in  the  following: 

V  P+  f  =0  VXeB0  and  Vf  e  [ 0,tf\  (45) 

* 

Based  on  the  definition  of  P  from  Equation  (6),  P  is  calculated  as 

P  =  det  F[tr(F  F  1  )T+  T  -T(F  F  1  )T]  F“T  (46) 


In  the  above  equation,  T  is  given  by  Equation  (43).  Suppose  that  the  gravity  is  the  only  body 
force  acting  on  the  body.  Let  po  be  the  constant  density  of  the  body  in  the  reference  configuration 
and  let  g  be  the  gravity  constant.  The  body  force  density  in  the  current  configuration  is  given  by 


b  =  pg  = 


detFs 
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Hence,  f  =  pog,  and  therefore. 


* 

f  =0 


(47) 


A  variational  form  of  the  kinematic  sensitivity  problem  can  now  be  posed  as  the  calculation  of 
x(X,  t)  such  that 


(V-P)-  u(X)dF  =  0 


(48) 


where  u(X)  is  a  kinematically  admissible  sensitivity  displacement  field  expressed  over  the  reference 
configuration.  With  integration  by  parts,  Equation  (48)  results  in 


*  30 


Pm  •  udH0 


(49) 


'0BO 


where  m  is  the  unit  outward  normal  at  X  6  <3B0.  The  reference  configuration  is  not  affected  by 

the  change  in  the  die  surface,  therefore  m  =0.  Transforming  the  term  in  the  right-hand  side  of 

Equation  (49),  to  the  current  configuration  (where  n  is  the  corresponding  unit  normal  vector)  and 

* 

substituting  the  expression  for  P  from  Equation  (46),  one  can  obtain. 


Ab0 


Pm  •  udH0  =  /  (tr(FF_1)T+  T  — T(FF“‘)T}n  •  udA„+\ 


(50) 


AB„ 


The  normal  vector  n  in  the  current  configuration  can  be  expressed  in  terms  of  the  deformation 
gradient  field  and  the  normal  vector  m  as  follows: 


n  =  F~Tm/||F-Tm| 


(51) 


Applying  the  sensitivity  operator  to  the  above  equation  and  after  some  simplifications,  one  can 
show  that 


n  =  {(F  F-1)  •  n  ®  n}n  -  (F  F_1)Tn 


(52) 


Using  the  above  equation,  the  sensitivity  of  the  traction  vector  t  =  Tn  (per  unit  current  area)  can 
be  written  as  follows: 


t  =  {(F  F-1)  •  n  ®  n}  t  +  Tn  —  T(F  F-1)  n 


(53) 


Equations  (49),  (50)  and  (53)  finally  lead  to  the  following  weak  sensitivity  problem:  Calculate 
x(X,  t)  such  that 


*  3u 


'0B„+i 


{(FF  ‘)  •  (I  -  n®  n)}t  •  udH„+i  + 


t  •  udA 


H+l 


(54) 


'0B„. 


for  every  admissible  test  function  u. 

In  the  direct  contact  analysis  presented  in  Sections  2.2  and  2.3,  the  contact  traction  k  per  unit 
reference  area  I  was  introduced.  The  contact  traction  vectors  t  and  k  are  related  as  follows: 


t  =  k||FrTn||/Jr 


(55) 
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where  Fr  =  FF„  is  the  relative  deformation  gradient  and  Jr  =  det  Fr.  Using  the  fact  that  |Fr  n||/./r  = 


d  -4  n 
d  An^-\ 


,  the  first  term  in  the  right-hand  side  of  Equation  (54)  can  be  expressed  as 


/  {(FF  1 )  •  (I  -  n  0  n)}  t  •  u &An+\  =  /  {(F  F  *)  •  (I  -  n  0  n)}X  •  udAn  (56) 

Jd B„+i  </eB„ 

The  sensitivity  of  the  contact  traction  vector  t  can  be  computed  as 

+  (57, 

Jr  Jr  Jr 

Equation  (57)  can  be  used  to  modify  the  right-hand  side  of  Equation  (54)  as  follows: 

[  i  -udAn+l=  [  +  U(L4,i+1  (58) 

iSB„+1  7oB„+i  \  Jr  Jr  Jr  J 

The  first  term  in  the  right-hand  side  of  the  above  equation  simplifies  as 

f  ^11  Ft  nil  /*  ^ 

/  X  y - ud^„+i=  /  X  -ud^„  (59) 

7SB„+1  Jr  J  SB,, 

* 

The  computation  of  the  sensitivity  of  the  contact  traction  vector  X  is  derived  in  the  next  section. 
The  following  expression  can  be  shown  to  hold: 


|FrTnj|  =  [(FrFr  )  •  (n  0n)  +  (Fr  F* )  •  (n®  n)]- 


1 


Fjnll 


where  the  sensitivity  of  the  relative  deformation  gradient  is  expressed  as 


Fr  =(F  — Fr  F„)F“ 


-l 


The  sensitivity  of  Jr  is  computed  as 

* 


detFr  =  Jrtr(F  F  x)  —  Jr  tr(Fr  F„  F  *) 


-r 


(60) 


(61) 


(62) 


Equations  (60)  and  (62)  can  be  used  to  evaluate  the  second  and  third  terms  in  the  right-hand  side 
of  Equation  (58).  One  can  thus  pose  the  weak  kinematic  sensitivity  problem  as  follows: 

Calculate  x(X,  t)  such  that 


£P  'dXdV°  =  . 


/  X  -ud^„+  /  tr[Fr  F„  F  *]X-ud^„ 

SB„  JdB„ 


[ 

FrF„  F-'FT 

l|F>||2 

(n  0  n)X  •  u  d An 


(63) 


In  the  above  equation,  the  integrations  over  the  boundary  c?B„  are  non-zero  only  on  its  subset 
r  where  contact  occurs. 
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Remark  2.  In  using  a  finite  element  method  to  solve  the  above  kinematic  sensitivity  problem 

for  x,  one  can  notice  that  the  term  on  the  left-hand  side  of  Equation  (63)  contributes  to  both 
the  stiffness  matrix  and  the  force  vector.  Considering  the  right-hand  side  of  the  equation  above, 
we  remark  that  its  first  term  contributes  to  both  the  stiffness  matrix  and  the  force  vector  (see 
Section  3.3),  whereas  its  second  and  third  terms  contribute  only  to  the  force  vector.  In  the  present 
continuum  sensitivity  analysis,  the  stiffness  matrix  for  the  sensitivity  problem  is  not  the  same  as 
the  linearized  tangent  stiffness  used  in  the  direct  kinematic  analysis,  even  though  the  relation  of 
*  * 

dP  with  dF  and  P  with  F,  respectively,  look  identical.  This  becomes  apparent  after  noticing  that 
*  * 

P  depends  on  T  (Equation  (46))  and  through  the  solution  of  the  constitutive  sensitivity  problem 
* 

T  depends  upon  the  history  of  the  sensitivities  of  the  material  state  and  Fe  (Equations  (42)  and 
(43)).  Such  dependencies  are  characteristic  of  the  sensitivity  problem  and  are  of  no  relevance  to 
the  direct  analysis. 

3.2.3.  Initial  conditions  and  non-contact  boundary  sensitivity  conditions.  The  following  set  of 
initial  conditions  are  used  in  the  direct  analysis: 

T(X,  0)  =  0,  s(X,0)  =  s0  (64) 

The  corresponding  initial  conditions  for  the  sensitivity  problem  are  the  following: 

T(X,0)  =  0,  *(X,0)  =  0  (65) 

For  calculating  sensitivities  with  respect  to  the  initial  state  or  stresses,  one  should  appropriately 
modify  these  equations. 

The  boundary  conditions  for  the  die  sensitivity  problem  are  derived  from  the  corresponding 
boundary  conditions  of  the  direct  analysis.  Some  typical  non-contact  related  cases  are  examined 
below.  The  boundary  traction  is  zero  on  a  free  boundary.  The  corresponding  boundary  condition 
for  the  sensitivity  problem  is  written  as  follows: 

t  =0  (66) 

On  the  part  of  the  boundary  with  prescribed  displacement  or  velocity,  e.g. 

x(X,0  =  x(X,0  (67) 

where  x(X,  t)  is  a  known  function,  the  corresponding  kinematic  sensitivity  condition  is 

x(X,t)  =  0  (68) 

The  equation  above  should  be  modified  if  the  objective  is  to  calculate  sensitivity  fields  with  respect 
to  design  variables  related  to  x(X,  t)  (for  example  when  the  ram  speed  is  the  design  variable). 

3.3.  The  contact  sensitivity  problem 

As  a  result  of  the  non-smooth  nature  of  the  contact/friction  conditions,  the  calculation  of  the  sensi¬ 
tivity  of  contact  traction  components  is  a  non-trivial  task.  Non-smooth  contact/frictional  behaviour 
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result  from  the  variation  of  the  normal  traction  component  at  impending  contact  and  from  the  vari¬ 
ation  of  tangential  traction  component  at  the  stick-slip  transition.  To  allow  for  the  differentiability 
of  the  direct  contact/frictional  conditions  (see  Section  2.3),  the  following  regularizing  assumptions 
are  made: 

(1)  Contact  sensitivity  assumption:  A  material  particle  that  lies  in  the  admissible  (inadmissible, 
respectively)  region  at  time  t  for  a  particular  deformation  problem  also  lies  in  the  admissible 
(inadmissible,  respectively)  region  for  a  nearby  deformation  problem  (e.g.  obtained  after  an 
infinitesimal  perturbation  of  the  reference  die  or  for  a  nearby  set  of  process  parameters)  at 
that  time  t. 

(2)  Friction  sensitivity  assumption:  If  the  frictional  conditions  at  a  particular  point  in  contact  for 
a  reference  deformation  problem  are  those  of  sticking  (sliding,  respectively),  the  frictional 
conditions  of  this  point  at  the  same  time  for  a  nearby  deformation  problem  will  be  that  of 
sticking  (sliding,  respectively).  Thus  transition  from  a  stick  to  a  slip  condition  does  not  occur 
at  a  material  point  as  a  result  of  infinitesimal  changes  in  the  die  shape  or  other  process 
parameters. 

Recall  that  in  the  direct  contact  problem,  in  addition  to  contact  tractions,  the  regions  of  contact 
and  the  regions  of  slip/stick  within  the  contact  boundary  are  calculated.  Flowever,  in  the  sensi¬ 
tivity  analysis  with  the  regularization  assumptions  introduced  above,  one  needs  to  calculate  only 
the  sensitivities  of  the  traction  components,  with  the  contact  boundary  and  regions  of  stick/slip 
identified  by  the  direct  analysis. 

Let  us  first  introduce  a  proper  notation  for  the  sensitivity  of  a  function  /  =  /(x).  In  the  particular 
case  of  a  die  design  problem,  it  is  assumed  that  a  variation  in  the  die  surface  would  not  only 

result  in  a  variation  of  the  position  x,  but  also  may  result  in  a  variation  of  the  function  /  itself. 

* 

The  sensitivity  of  the  function  /  denoted  by  [/(x)]  generally  consists  of  two  parts,  the  first 
* 

part,  / ,  denotes  the  contribution  due  to  changes  in  the  function  /  with  x  constant  and  the  second 
part,  V/x,  denotes  the  contribution  due  to  changes  in  the  variable  x  with  /  constant,  i.e. 

L/fel  =  /  +  V/x  (69) 


A  calculation  of  the  sensitivity  X  of  the  contact  traction  is  now  developed.  Figure  3  shows  the 
schematic  used  to  define  the  contact  traction  sensitivities.  The  analysis  presented  here  is  for  die 

design  problems  and  the  change  in  the  die  shape  is  represented  by  y. 

* 

In  the  present  work,  the  sensitivity  traction  field  X  is  defined  by  differentiating  with  respect 
to  the  die  surface  the  continuum  contact  constraints  and  then  implementing  these  sensitivity  con¬ 
tact  constraints  consistently  with  the  analysis  used  in  the  direct  contact  problem.  The  variables 

(XnAn)  are  assumed  to  be  known  as  calculated  from  the  time  integration  of  the  contact  sensitiv- 
ity  problem  from  the  previous  time-step  (time  t„).  The  corresponding  unknown  variables  at  time 
*  * 

4+i  are  denoted  as  (X,c ),  where  for  simplicity  of  notation,  the  subscript  (n  +  1)  is  omitted  from 
variables  that  refer  to  time  tn+  \ .  The  solution  of  the  direct  deformation  contact  problem  gives  the 
contact  tractions  X  at  time  4+1  and  identifies  regions  of  sticking  and  sliding. 
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Perturbed 

Die 

•i* 

y  =  y  +  y 


Figure  3.  A  schematic  of  the  contact  sensitivity  problem  for  infinitesimal  changes  in  the  die  shape. 


Consideration  of  the  strong  form  of  the  normal  contact  constraints  in  the  direct  deformation 
problem  yields  for  points  in  contact  gf(x„+i)  =  0  and  > 0.  Thus,  the  strong  form  of  the  normal 
contact  sensitivity  constraint  is  obtained  as 

_ 5j« 

gf(x„+i )  =  0  and  An  €  3? 

* 

This  is  a  direct  result  of  the  contact  sensitivity  assumption  introduced  earlier.  The  sensitivity  2n 
is  treated  as  the  Lagrange  multiplier  that  enforces  the  continuum  form  of  the  normal  contact 
sensitivity  constraint.  The  following  penalty  form  is  used  to  enforce  this  constraint: 

*  *  — — — 

7n  =  An„  +£n0(x„+i)  (70) 

* 

where  An,,  is  the  sensitivity  of  the  normal  traction  component  at  the  beginning  of  the  time  integra¬ 
tion  step.  The  penalty  parameter  eN  introduced  above  is  generally  not  the  same  as  the  corresponding 
parameter  used  in  the  direct  contact  sub-problem  (see  Remark  3). 

Consideration  of  the  tangential  contact  constraints  yields  that  for  points  in  sticking  contact, 
|  =  0.  The  corresponding  sensitivity  constraint  (based  on  the  friction  sensitivity  assumption)  takes 
the  form 

*  * 
l  =  1=0 

* 

As  in  the  case  of  the  normal  to  the  die  sensitivity  constraints,  At  is  treated  as  the  Lagrange 
multiplier  that  enforces  the  tangential  to  the  die  sensitivity  constraints  for  the  case  of  sticking 
contact.  The  following  penalty  formulation  of  the  sticking  sensitivity  constraint  is  introduced: 

*  1  * 

l  =  -IT  (71) 

fiT 
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integration  of  this  equation  results  in  the  following: 

^  ^  *  * 

At  =  At„  +£t(I  -  l„)  (72) 

* 

where  /'t„  is  the  sensitivity  of  the  tangential  traction  component  at  the  beginning  of  the  time 
integration  step.  The  parameter  et  can  be  selected  independently  of  the  corresponding  parameter 
in  the  direct  contact  sub-problem. 

Remark  3.  Equations  (70)  and  (72)  are  derived  in  a  continuum  contact  sensitivity  setting.  An 
alternate  formulation  can  be  considered  where  the  design  differentiation  of  the  corresponding  dis¬ 
crete  equations  used  in  the  augmented  Lagrangian  analysis  in  the  direct  contact  sub-problem  is 

carried  out  (i.e.  of  An  =  (An  +  £n0(xJ,+i))  and  /U[.ial  =  AT„  +  £j(  11+ i  —  f„)  +  AAlj\  respectively). 

However,  this  interpretation  requires  that  the  penalty  parameters  en  and  et  are  identical  in  both 
the  direct  and  sensitivity  contact  sub-problems.  Such  a  restriction  is  unnecessary  as  it  implies  that 
the  magnitude  of  the  penalty  parameters  en  and  et  in  the  sensitivity  contact  analysis  is  limited 
by  corresponding  values  used  in  the  direct  contact  problem.  Indeed,  recall  that  an  augmented 
Lagrangian  formulation  instead  of  a  penalty  formulation  is  used  in  the  direct  contact  algorithm 
because  using  large  penalties  to  enforce  contact  and  frictional  constraints  leads  to  ill-conditioning 
and  convergence  problems  of  the  Newton  algorithm  that  solves  the  non-linear  direct  deformation 
boundary  value  problem. 

Remark  4.  Solution  of  the  direct  deformation  problem  with  moderate  penalties  within  a  time 
increment  usually  requires  augmentations  to  accurately  enforce  the  contact/frictional  constraints. 
Derivation  of  the  contact  sensitivity  sub-problem  by  differentiation  with  respect  to  the  die  of 
the  discrete  direct  contact  sub-problem  (see  Remark  3)  will  thus  lead  to  a  contact  sensitivity 
algorithm  that  requires  a  number  of  iterations  (augmentations).  The  sensitivity  deformation  problem 
is  however  a  linear  problem  and  augmentations  are  preferably  avoided  with  a  linear  continuum 
sensitivity  analysis. 

As  a  result  of  the  contact  and  frictional  sensitivity  assumptions.  Equations  (70)  and  (72)  can 
be  used  here  with  high  penalty  parameters  and  no  augmentations  are  necessary  in  the  contact 
sensitivity  analysis.  For  simplicity  of  notation,  the  contact  sensitivity  penalty  parameters  are  denoted 
by  en  and  et,  even  though  they  are  not  the  same  as  the  corresponding  penalty  parameters  in  the 
direct  contact  analysis. 

* 

In  the  case  of  slip.  At  is  well  defined  by  the  Coulomb  friction  law  as  follows: 

_ * _ 

At  =  (73) 

To  simplify  the  above  equation,  the  following  relation  is  used: 

*  *  *  - 

[H(y)J  =  [y^I  =  y,{  +y,«  l  (74) 

* 

where  y  =  represents  the  contribution  to  the  sensitivity  due  to  changes  in  the  shape  y  of  the  die. 
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Using  the  above  equation.  Equation  (73)  (applicable  to  a  slip  condition)  simplifies  as  follows: 

(75) 


*T  =  ^at-at 

An  \  y,{  •  y,c 


y,i  •  ?,{ 


_AT|Li^2|{ 


y,{  •  y,c 


One  can  thus  complete  the  contact  sensitivity  description  by  writing  the  following  equation  for  X: 

^  *L-  ^ 


X  =  Xn  ~  Xt ,  Xn  —  ^Nv(y),  It  —  ^Tti(y) 


(76) 


where  Xn  and  Xt  are  the  sensitivities  of  the  normal  and  tangential  contact  traction  vectors, 
respectively.  Using  Equations  (70),  (74),  (76b)  and  v(y)  =  ti(y)  x  e2/||Ti(  y)||,  the  sensitivity 
of  the  normal  contact  traction  vector  can  finally  be  expressed  as  follows: 


An  —  dC,l+  \  C,2  +  £3 


(77) 


where  the  vectors  C;  are  defined  by 

(y,i  x  e2) 


Cl  =  CN- 


|y,fl 


C2  = 


An 

I  y,s\ 


(y^xe2)-(^^)(y;{xe2) 


(78) 


-  *  (y,{  x  e2)  ,  An 

S3  ii  —  M  i 


|y,d 


y,ii 


(y,f  xe2)  — 


* 

y,g  •  y,{ 
y.c  •  y.« 


(fj;  x  e2) 


Similarly,  using  Equations  (72)  (for  sticking  friction)  or  (75)  (for  sliding  friction),  (74)  and  (76c), 

* 

the  sensitivity  of  the  tangential  contact  traction  vector  Xt  can  be  expressed  as  follows: 

*  * 


Xt  —  9  3 1  +  £,  32  +  $3 

If  the  material  point  is  in  sliding  contact,  then  the  vectors  3,  take  the  form: 


(79) 


£n  . 

*  At 

An 

At^  -  ( 

'y,f-y,« 

v  irlt 

* 

An"  Xt  - 

(  y ,4  •  h 

,  at 

An 

{ y,{*y,{ 

XT  +  At  y,£ 


(80) 
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On  the  other  hand,  if  the  friction  conditions  are  that  of  stick,  3,  take  the  form 

31  =  0 

32  =  £t  +  At  (81 ) 

Jj 4  ^ 

a3  =  At„  y>{  -  &r  y>{  +  At  y,c 

If  the  contact  conditions  are  frictionless,  then  3,  =  0. 

Substitution  of  Equations  (77)  and  (79)  in  Equation  (76a)  yields  a  concise  expression  for  the 
sensitivity  of  the  contact  tractions  as 

*  *  * 

k  =  0  71+ f  72  +  73  (82) 

.  -  * 
where  y,-  =  C,  —  3,.  The  sensitivity  c  of  the  amount  of  inelastic  slip  is  related  to  x  (see  Figure  3). 

-  * 

The  linear  relationship  between  c  and  x  is  now  developed.  Consider  the  sensitivity  of  the  relation 
( y  —  x)  •  ti(  y)  =  0  which  expresses  the  fact  that  the  closest  point  y  on  the  die  to  a  point  x  is  the 
projection  of  the  point  x  on  to  the  die.  The  sensitivity  of  the  above  relation  is  given  by 

*  * 

(y(l)-x).[y^(O]+([y(O]-x).y;{(|)  =  0  (83) 

The  above  expression  can  be  simplified  as  follows: 

I  =  a  •  x  +b  (84) 


where 

a-  Ti(y> 

l|ih(y)||2[i  +0K(y)] 

_  *  _  * 
~[gv(y)-  y,g  +ti(y)-  y] 

IKi(y)ll2[l  +3»e(y)] 


(85) 


where  k(j)  =  v(y)  •  y ,^(|)/||Ti(y)||2  represents  the  curvature  of  the  die  at  y. 

*  * 

The  linear  relationship  between  the  sensitivity  of  the  gap  function  g(x)  and  x  can  be  obtained 
by  the  design  differentiation  of  Equation  (19)  as  follows: 

9  =v(y)-(|  -  x)  (86) 

Substitution  of  Equations  (84)  and  (85)  in  Equation  (82)  results  in  the  final  expression  for  the 
sensitivity  of  the  contact  traction: 


l  =  [(72  0  a)-  (71  0  v(y))]  x  +[(v(y)  •  y)7i  +  by2  +  73] 


(87) 


From  Equations  (63)  and  (87),  one  can  conclude  that  the  sensitivity  of  the  contact  tractions 
contributes  to  the  force  vector  as  well  as  the  stiffness  matrix  in  the  weak  form  of  the  kinematic 
sensitivity  problem  (see  the  first  term  on  the  right-hand  side  of  Equation  (63)). 
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Figure  4.  Solution  procedure  used  in  the  contact  sensitivity  sub-problem.  The  penalty  parameters  en 
and  £t  used  in  this  analysis  are  much  larger  than  those  used  in  the  direct  contact  sub-problem  in  which 

augmentations  are  usually  necessary. 

The  initial  conditions  for  the  contact  sensitivity  algorithm  at  the  beginning  of  the  deformation 

'I'  — 

at  time  to  are  (ko,  c0)  which  are  given  as  follows: 

*  * 

ko  =0,  l0  =0  (88) 

A  summary  of  the  time  integration  of  the  contact  sensitivity  sub-problem  follows  (Figure  4): 
Solve  for  the  sensitivity  of  the  deformation  x  at  time  tn+\ 

5(x,  u)  =  S(x,  u)  +  Sc(x,  u)  =  0 

Note:  The  above  equation  represents  the  linear  sensitivity  weak  form,  the  solution  of  which  re¬ 
quires  input  from  the  kinematic,  constitutive  and  contact  sensitivity  sub-problems.  The  solution  of 

this  linear  equation  yields  the  sensitivity  of  the  body  configuration  x  corresponding  to  the  body 
configuration  B„+i. 

The  following  values  for  the  sensitivity  of  contact  tractions  are  used  in  the  solution  for  ma¬ 
terial  points  in  contact  (based  on  the  computed  body  configuration  Bh+i  and  the  history  of  the 
sensitivity): 

Sensitivity  of  normal  traction: 

* 

Use  >,k  from  Equation  (77) 

Sensitivity  of  tangential  traction: 

IF  (xeystick)  THEN 
* 

Use  It  from  Equations  (79)  and  (81) 

ELSE  (sliding  contact) 

* 

Use  It  from  Equations  (79)  and  (80) 
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Remark  5.  Implementation  of  the  contact  sensitivity  algorithm  obtained  by  differentiation  of  the 
discrete  equations  used  in  the  direct  contact  sub-problem  was  also  developed  as  part  of  this  in¬ 
vestigation.  Such  a  contact  sensitivity  analysis  uses  the  same  penalty  parameters  sN  and  ex  as  the 
direct  analysis  and  requires  augmentations  within  a  time  increment.  It  has  been  determined  that 
the  method  presented  in  this  paper  provides  the  same  accuracy  for  the  computed  sensitivities  of 
the  contact  tractions  as  this  alternative  method,  although  with  a  higher  computational  efficiency. 


4.  NUMERICAL  RESULTS 

A  comprehensive  set  of  numerical  examples  in  the  design  of  metal  forming  processes  is  considered 
here  to  validate  the  proposed  computational  framework  and  numerical  algorithm.  The  CPU  require¬ 
ments  of  the  present  object  oriented  implementation  are  shown  in  a  representative  example.  The 
computations  were  performed  on  an  IBM  RS-6000  workstation  at  the  Cornell  Theory  Center.  In 
all  reported  simulations  three-noded  cross  triangular  elements  were  used.  Unless  otherwise  stated, 
the  penalty  parameters  for  the  sensitivity  contact  problem  are  selected  as  eN  =  108  and  i;T  =  1 04. 
The  workpiece  material  chosen  in  all  examples  is  A1  1100-0  at  673  K.  The  constitutive  model  for 
this  material  is  taken  from  Brown  et  al.  [27]  and  is  summarized  in  Appendix  B. 

4.1.  Accuracy  investigations 

4.1.1.  Sensitivity  calculations  in  axisymmetric  extrusion  ( Example  1).  An  axisymmetric  extrusion 
of  a  cylindrical  workpiece  is  considered  through  a  curved  die.  The  initial  radius  of  the  workpiece 
is  1  cm  and  the  initial  height  is  2  cm.  A  friction  coefficient  of  0.01  is  assumed  at  the  die-workpiece 
interface. 

A  degree  six  («  =  6)  Bezier  curve  is  used  to  represent  the  die  shape: 


n+\ 

r  =  E  Q/K*) 

;= l 

z  =  a  (in  cm) 


(89) 

(90) 


Here  a  e  [0, 1],  and  C„  i  =  1, . . . ,(«+ 1 ),  are  algebraic  control  parameters  (points).  The  Bernstein 
functions  /}(a),  /=!,...,(«  +  1),  are  defined  as  follows: 


/i(a) 


n\ 


j- 1 


(i  -  !)!(«  -  i  +  1)! 


(1 


a) 


n—i+ 1 


(91) 


In  order  to  obtain  the  same  reduction  for  different  die  design  parameters,  the  parameters  C, 
and  C7  are  selected  such  that  C,  =  1.0  cm  and  C7  =  0.866  cm.  Such  a  selection  results  in  a  fixed 
cross-sectional  area  reduction  of  25  per  cent.  In  addition,  to  allow  the  workpiece  to  enter  and  leave 
the  die  zone  smoothly,  the  die  curve  at  the  entrance  and  exit  is  restricted  to  be  perpendicular  to 
the  radial  axis.  This  is  accomplished  by  selecting  C2  =  Cj  and  C(,  =  Cj. 

With  the  above  selection  of  the  parameters  C,,  there  are  three  remaining  die  design  parameters 
(/li, ^2,163)  =  (U3, C4,  U5 ).  The  initial  values  are  arbitrary,  and  the  reference  (initial)  die  is  selected 
such  that 


II,  =  (3Cj  +  C7)/4,  p2  =  ( C j  +  C7)/2,  /1 3  =  (Cj  +  3C7)/4 
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Figure  5.  The  initial  and  steady-state  (7=  170  s)  Figure  6.  Steady-state  distribution  of  the  state 
grids  shown  along  with  the  initial  guess  of  variable  s  (in  MPa)  (Example  1). 

the  die  shape  (Example  1). 


Table  1.  Simulation  parameters  for  the  frictional 
extrusion  of  a  1100-A1  billet  (Example  1). 


Parameter 

Value 

Energy  error  norm 

1.0E-07 

Displacement  L2  error  norm 

1.0E-07 

Normal  penalty  for  contact 

1.0E  +  06 

Tolerance  for  gap 

1.0E-06 

Tangent  penalty  for  contact 

1.0E  +  03 

Tolerance  for  friction  condition 

1.0E-06 

The  initial  die  shape  (length=lcm)  is  shown  in  Figure  5.  A  finite  element  mesh  of  6  x  18 
cross-triangular  elements  is  considered  (Figure  5).  The  extrusion  velocity  is  0.01  cm/s.  A  fixed 
time  step  of  At  =  0.25  s  is  used.  The  parameters  used  in  the  simulation  are  given  in  Table  1.  The 
time  at  which  the  process  reached  steady-state  conditions  at  the  exit  was  computed  as  t  =  1 70  s. 
The  grid  at  this  time  is  also  shown  in  Figure  5. 

Figure  6  shows  the  steady-state  contour  plot  of  the  state  variable  s.  In  more  detail,  the  state 
variable  distribution  s  along  the  die  exit  is  shown  in  Figure  7  demonstrating  that  steady-state 
conditions  are  achieved  at  1=  170  s.  Figures  8-10  show  the  steady-state  contours  of  the  sensitivity 

5  with  respect  to  the  design  parameters  (/?i ,  /?2,  /?3 )  using  the  DDM  (direct  differentiation)  and 
the  FDM  (finite  difference)  methods.  The  FDM  results  are  obtained  using  the  results  of  the  direct 
analysis  and  a  forward  difference  approximation  for  a  perturbation  of  each  of  the  design  variables 
by  10_Jcm.  Figure  11  shows  a  comparison  of  the  sensitivity  s  with  respect  to  each  of  the  design 
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r(cm) 


Figure  7.  The  distribution  of  the  state  variable  s  along  the  exit  at  various  process  times  (Example  1). 


Figure  8.  Contours  of  the  steady-state  sensitivity  of  the  state  variable  s  (in  MPa)  with  respect  to  (h  computed 

using  the  DDM  and  FDM  methods  (Example  1). 


parameters  along  the  die  exit  using  the  DDM  and  the  FDM  at  steady  state.  The  results  match 
quite  well. 

4.1.2.  Evaluation  of  the  accuracy  of  the  computed  sensitivities  of  contact  tractions  in  an  open 
die  forging  process  ( Example  2).  In  this  example,  the  accuracy  of  the  computed  sensitivities  of 
contact  tractions  with  respect  to  the  speed  of  a  flat  die  in  an  open  die  forging  process  is  examined, 
i.e.  the  forging  rate  V  is  here  taken  as  the  process  (design)  parameter.  Variations  in  the  forging 
rate  affect  the  dissipated  energy  during  upsetting  and  so  it  is  important  to  be  able  to  model  the 
sensitivity  fields  as  a  result  of  this  variation. 

The  initial  cylindrical  workpiece  of  radius  1.00mm  and  height  3.00mm  is  discretized  using 
576  triangular  finite  elements.  As  a  result  of  the  symmetry  of  the  problem,  only  a  quarter  of  the 
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Figure  9.  Contours  of  the  steady-state  sensitivity  of  the  state  variable  s  (in  MPa)  with  respect  to  fh  computed 

using  the  DDM  and  FDM  methods  (Example  1). 


Sensitivity  of  s 
0.008 
0.0045 
0.001 
-0.0025 
-0.006 
-0.0095 
-0.013 
-0.0165 
-0.02 


Figure  10.  Contours  of  the  steady-state  sensitivity  of  the  state  variable  s  (in  MPa  )  with  respect  to  (h  computed 

using  the  DDM  and  FDM  methods  (Example  1). 


work-piece  is  modelled  (Figure  12).  The  specified  forging  velocity  for  the  reference  problem  is 
V  =  0.01  mm/s  and  the  workpiece  is  allowed  to  deform  to  80  per  cent  of  its  initial  height  at  which 
time  the  sensitivity  fields  are  monitored.  The  various  simulation  parameters  are  given  in  Table  11. 
The  FDM  values  are  computed  corresponding  to  a  variation  of  SV  =  10_4mm/s. 

The  accuracy  of  the  contact  sensitivity  algorithm  introduced  in  this  paper  (see  Figure  4)  is 
investigated  here  with  the  penalty  parameters  in  the  sensitivity  algorithm  selected  as  eN  =  1.02s +08 
and  eT  =  1 .02s  +  06.  In  addition,  an  alternative  simulation  is  also  performed  with  smaller  contact 
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Figure  11.  Comparison  of  steady-state  values  of  s  at  the  exit  corresponding  to  /?i ,  [>2  and  [h ,  respectively, 
computed  using  the  DDM  and  FDM  methods  (Example  1). 


'  j  t3ic  !  v 


Asis  of  cylinder 


Eire  ■■'jrljL'.' 


T 

F  i | ■  1 . 111.1l  in.” 


Figure  12.  Axisymmetric  upset  forging  for  parameter  sensitivity  analysis 
with  respect  to  the  die  speed  (Example  2). 


penalty  parameters.  As  mentioned  in  Remarks  4  and  5,  augmentations  are  required  within  a  time 
step  for  such  a  choice.  The  penalty  parameters  for  this  simulation  are  selected  to  be  the  same 
as  the  corresponding  parameters  in  the  direct  contact  problem.  Such  an  alternative  simulation  can 
be  thought  of  as  the  regularized  derivative  of  the  corresponding  discrete  contact  algorithm  of  the 
direct  analysis. 
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Table  II.  Simulation  parameters  for  the  fric¬ 
tional  open  die  forging  of  an  axisymmetric 
1100-A1  billet  (Example  2). 


Parameter 

Value 

Energy  error  norm 

1.0E  - 

08 

Displacement  7,2  error  norm 

1.0E  - 

08 

Normal  penalty  for  contact 

1.0E  + 

07 

Tangent  penalty  for  contact 

1.0E  + 

05 

Tolerance  for  gap 

1.0E  - 

08 

Tolerance  for  friction 

1.0E  - 

06 

50  - 1 — i — i — | — i — i - 1 — | — i — i - 

40  1  *  Normal  traction 

;  6  Tangential  traction 

_  30  - 
cri 


0.0  0.2  0.4  0.6  0.8  1.0 


Radius  at  the  die-workpiece  interface  {  mm  ) 


Figure  13.  Distribution  of  the  nonnal  and  tangential  contact  traction  components  at  the  final  state  for  the 
axisymmetric  frictional  upsetting  process  with  a  reference  die  velocity  of  0.01  mm/s  (Example  2). 


In  order  to  validate  the  parameter  sensitivity  problem  for  processes  involving  frictional  contact, 
a  moderate  friction  coefficient  of  0.4  is  chosen  to  allow  for  variable  regions  of  stick  and  slip  at 
the  contact  interface.  Figure  13  depicts  the  variation  of  the  normal  and  tangential  tractions  at  the 
die-workpiece  interface  using  the  reference  die  velocity.  Figure  14  shows  the  variation  of  dX^/dV 
and  dXj/dV  along  the  radius  in  the  deformed  configuration  at  the  die-workpiece  contact  interface. 
In  general  the  FDM  and  DDM  results  (using  the  present  contact  sensitivity  algorithm  and  the 
alternative  contact  sensitivity  algorithm  mentioned  above)  agree  very  well.  It  is  observed  that  the 
DDM  method  is  able  to  accurately  predict  the  tangential  traction  derivatives  in  regions  where  there 
is  a  transition  from  stick  to  slip  or  vice  versa.  This  is  evident  at  the  sharp  peak  at  radius  r  ~  0.7 
in  Figure  14  where  there  is  a  transition  from  sticking  friction  (inner  radii)  to  sliding  friction  (outer 
radii).  The  non-smooth  behaviour  at  this  location  is  the  result  of  the  Coulomb  law  used  to  model 
friction.  Small  oscillations  in  the  sensitivities  (observed  both  in  the  FDM  and  DDM  methods)  are 
due  to  the  instability  in  the  finite  element  implementation  of  the  direct  deformation  problem  as  a 
result  of  the  reduced  numerical  quadrature  used  in  the  treatment  of  near-incompressibility  and  is 
currently  being  investigated. 

The  computational  savings  using  the  present  DDM  method  are  substantial  for  problems  where 
the  contact  conditions  are  complex  and  is  therefore  the  recommended  method  that  is  used  in  the 
remaining  of  this  paper. 
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Radius  at  the  die-workpiece  interface  {  mm  ) 


Figure  14.  Comparison  of  the  design  derivatives  ak^/tlV  and  dkj/dV  computed  with  the  DDM  and  FDM 
methods  for  the  axisymmetric  frictional  upsetting  process  (Example  2). 


Die 


workpiece 


Figure  15.  Die  design  for  uniform  distribution  of  the  state  variable  s  at  the  exit  location 
in  axisymmetric  extrusion  (Example  3). 

4.2.  Using  sensitivity  calculations  for  die  design  problems 

4.2.1.  Die  design  for  uniform  material  state  in  an  extruded  product  ( Example  3).  In  this  exam¬ 
ple,  the  die  parameters  in  the  Bezier  representation  given  in  Section  4.1.1  are  designed  such  that 
the  distribution  of  s  at  the  die  exit  (see  Figure  15)  is  as  uniform  as  possible,  i.e.  one  is  interested 
in  calculating  the  die  parameters  such  that  the  following  function  is  minimized: 

/(M  =  £(*-*)2  (92) 

(=1 

where  (i=l,...,N)  are  N  (=20)  discrete  (equally  spaced)  state  variable  s  values  along  the 
radius  at  the  exit  (in  the  particular  example,  along  z=  1.0cm)  and  s=  is  the  mean 

value  of  the  distribution  of  s  along  the  exit. 

The  extrusion  conditions  are  identical  to  those  given  in  Example  1.  The  initial  mesh  is  shown 
in  Figure  5  with  6x18  cross  triangular  elements.  A  fixed  time  step  A t  =  0.25  s  is  considered  and 
the  steady-state  process  time  is  taken  as  1=  170  s. 

The  computed  sensitivity  fields  are  used  to  evaluate  the  gradient  of  the  objective  function.  Vari¬ 
ous  optimization  methods  have  been  tested  for  this  design  problem.  Figure  16  shows  the  objective 
function  values  as  well  as  the  norm  of  the  gradient  of  the  objective  function,  during  the  iterations 
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Figure  16.  The  objective  function  and  the  norm  of  gradient  of  the  objective  function,  respectively, 

versus  the  iteration  index  (Example  3). 


Figure  17.  Distribution  of  s  along  the  exit  using  the  SD  and  CG  methods,  respectively  (Example  3). 


using  the  steepest  descent  (SD)  and  conjugate  gradient  (CG)  methods.  Figure  17  shows  the  distri¬ 
bution  of  s  along  the  exit  at  the  final  time  during  various  optimization  iterations.  Convergence  is 
achieved  after  14  iterations.  The  computed  die  shapes  during  the  initial  and  following  optimization 
iterations  using  the  SD  and  CG  methods  are  shown  in  Figure  18.  Both  optimization  methods  result 
in  the  same  optimal  die. 

The  computing  statistics  for  the  optimization  of  the  die  shape  in  one  iteration  are  summarized 
in  Table  111.  In  the  FDM  method,  one  optimization  iteration  includes  4  direct  problems.  In  the 
DDM  method,  one  optimization  iteration  includes  1  direct  problem  and  3  sensitivity  sub-problems 
The  DDM  method  results  in  savings  of  about  64  per  cent  compared  to  the  FDM  method. 

4.2.2.  Die  design  for  minimum  required  extrusion  force  for  a  given  reduction  ratio  ( Example  4). 
The  same  Bezier  approximation  and  design  parameters  are  used  as  in  the  earlier  extrusion  examples. 
The  die  design  parameters  are  calculated  such  that  the  required  extrusion  force  is  minimized. 
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Figure  18.  The  computed  die  shapes  during  various  iterations  using  the  SD 
and  CG  methods,  respectively  (Example  3). 


Table  Ill.  CPU  time  (in  hours)  for  the  extrusion 
die  design  problem  (Example  3). 


Direct  problem  2.5 

Sensitivity  sub-problem  0.3 

1  optimization  iteration  with  FDM  10.0 

1  optimization  iteration  with  DDM  3.6 


The  required  force  is  expressed  as  follows: 

Fe=  [  l  ■  ez  dAn  (93) 

JBB„ 

* 

and  the  sensitivity  Fc  can  be  expressed  as 

Fe  =  l  -ezdAn  (94) 

J  3B„ 

where  the  integrations  over  the  boundary  dB„  are  non-zero  only  on  the  boundary  I  c  dB„  where 
contact  occurs. 

The  initial  design  parameters  are  selected  as  the  optimal  parameters  evaluated  in  Example  3,  i.e. 
pi  =  0.90829  cm,  fl2  —  0.88793  cm,  and  />3  =  0.88279  cm.  Due  to  the  finite  element  implementation 
of  the  contact  problem,  the  computed  extrusion  force  for  a  given  die  is  oscillatory  (Figure  19). 
Flowever,  when  the  stroke  is  greater  than  1.2  cm,  the  mean  extrusion  force  remains  constant.  The 
average  force  between  stroke  1.2  and  1.7  cm  is  used  here  as  the  objective  function  to  be  minimized. 
Figure  20  shows  the  sensitivity  of  the  extrusion  force  versus  stroke  with  respect  to  the  three  die 
design  parameters,  using  both  the  FDM  and  the  DDM  methods.  The  results  agree  quite  well.  In 
the  FDM  calculations,  the  increments  A/i;,  (=1,2,3,  were  taken  as  10_3cm. 

A  quasi-Newton  method  with  a  BFGS  update  of  the  Flessian  is  used  to  minimize  the  steady-state 
extrusion  force.  The  objective  function  is  approximated  at  each  optimization  step  as  a  quadratic 
function  along  the  descent  direction.  Three  different  step  sizes  are  selected  along  the  descent 
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Figure  19.  The  extrusion  force  versus  stroke  for  the  initial  die  design  (Example  4). 


Figure  20.  The  sensitivity  of  the  steady-state  extrusion  force  with  respect 
to  Pi,  P 2  and  (h,  respectively  (Example  4). 
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Figure  21.  The  extrusion  force  versus  stroke  and  the  various  die  shapes,  respectively,  computed  during 

the  optimization  process  (Example  4). 


direction  and  the  direct  problem  is  simulated  with  the  corresponding  design  variables  to  obtain 
the  values  of  the  objective  function.  The  three  computed  values  of  the  objective  function  define  a 
local  quadratic  approximation  that  can  be  used  to  find  the  minimum  in  the  descent  direction.  The 
extrusion  force-stroke  curve  and  the  corresponding  die  shapes  during  the  optimization  processes 
are  shown  in  Figure  21.  The  optimized  die  is  quite  flat  in  comparison  to  the  initial  die. 

4.2.3.  Die  design  in  axisymmetric  open-die  forging  to  obtain  a  product  of  desired  shape  ( Example 
5).  Consider  an  arbitrary  curved  die  with  which  an  axisymmetric  cylindrical  workpiece  is  forged 
for  a  given  reduction  resulting  in  a  product  of  some  final  shape.  In  this  example,  we  would  like  to 
use  an  optimization  scheme  with  the  present  sensitivity  analysis  to  recover  the  die  shape  starting 
from  the  shape  of  the  free  surface  of  the  final  product.  One  could  start  from  the  whole  shape  of 
the  final  product.  However,  since  the  top  surface  in  principle  defines  the  desired  die,  only  the  free 
surface  is  used  in  the  present  calculations.  Such  calculations  will  proceed  in  an  iterative  manner 
starting  with  a  flat  die  until  a  die  is  found  that  results  in  the  desired  shape  of  the  free  surface 
in  the  final  product.  A  single  stage  forging  process  cannot  always  produce  a  product  of  desired 
shape  and  generally  a  multi-stage  process  is  needed.  However,  based  on  the  direct  analysis  that 
produced  the  desired  shape  in  the  current  application,  it  is  known  that  a  die  does  exist  that  results 
in  a  product  of  the  desired  free  surface  shape  after  open-die  forging  of  the  workpiece. 

The  initial  radius  of  the  workpiece  is  1  cm  and  the  initial  height  3  cm.  A  friction  coefficient  of 
1.0  is  assumed  along  the  die-workpiece  interface.  A  degree  four  ( «  —  4 )  Bezier  curve  is  used  to 
represent  the  die  shape: 


n+ 1 

r=  1.5a  (in  cm),  z  =  )P  C,//(a) 

i=  1 

Here  ae[0, 1],  and  C„  i  =  1, ...,(«  +  1),  are  algebraic  parameters  with  f(a),  /  =  1, ...,(«  +  1), 
the  Bernstein  basis  functions  (see  Equation  (91)).  To  account  for  the  axisymmetry  of  the  problem, 
the  parameters  Ci  and  C2  are  selected  as  follows  C\  =  C2  (=  1.8  cm).  A  particular  die  is  defined 
with  (/i^ ,  ^2,  /?3 )  =  (C3,  C4,  C$ )  with  initial  values  of  (2.04, 2.13, 2.05)  cm.  This  die  shape,  the  initial 
workpiece  and  the  final  product  resulting  after  33.33  per  cent  height  reduction  (at  r  =  0)  are  shown 
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Figure  22.  The  initial  workpiece  and  selected  die.  After  33.33  per  cent  height  reduction  (at  r  =  0),  the  product 
shown  in  the  right  of  the  figure  results.  The  shape  of  the  free  surface  of  this  product  is  used  as  the  desired 

shape  in  the  optimization  analysis  (Example  5). 


Table  IV.  Simulation  parameters  for  the  design  of 
an  open  die  forging  process  (Example  5). 


Parameter 

Value 

Energy  error  norm 

1.0E  -  07 

Displacement  L  -x  error  norm 

1.0E  -  07 

Normal  penalty  for  contact 

1.0E  +  05 

Tangent  penalty  for  contact 

1.0E  +  04 

Tolerance  for  gap 

1.0E  -  05 

Tolerance  for  friction  condition 

1.0E  -  05 

in  Figure  22.  The  die  speed  is  taken  as  0.01  cm/s.  Due  to  axisymmetry,  only  1/4  of  the  workpiece 
cross-section  is  modelled  using  a  mesh  of  8  x  8  cross-triangular  finite  elements.  The  time  step 
was  selected  as  At  =  0.2  s.  The  remaining  parameters  used  in  this  direct  simulation  are  given  in 
Table  IV. 

In  the  present  design  problem,  one  is  given  the  initial  shape  of  the  workpiece  (Figure  22(a))  and 
is  interested  in  finding  the  die  shape  that  after  a  33.33  per  cent  height  reduction  (at  r  =  0)  results 
in  a  free  surface  of  the  final  product  as  that  given  in  Figure  22(b).  It  is  obvious  that  the  solution 
to  this  problem  corresponds  to  the  curved  die  defined  earlier.  The  CG  method  is  used  to  calculate 
the  descent  direction  during  iterations  and  a  line  search  method  is  used  to  calculate  the  step  size 
in  each  optimization  step.  Figure  23  shows  the  final  workpiece  and  the  computed  dies  during  the 
optimization  process.  Figure  24  shows  the  values  of  the  objective  function  as  well  as  the  norm  of 
the  gradient  of  the  objective  function  versus  the  number  of  iterations.  It  is  apparent  that  after  four 
iterations,  the  die  shape  has  converged  to  the  expected  solution  of  this  design  problem. 

Such  a  die  design  problem  is  essential  in  multi-stage  process  design  in  which  various  inter¬ 
mediate  preforms  have  to  be  designed.  Such  applications  will  be  reported  in  a  forthcoming 
publication  [34], 
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Figure  23.  The  final  workpiece  during  various  steps  in  the  optimization  process.  The  initial,  1st, 
2nd,  3rd  and  final  (4th)  iterations  are  shown  together  with  the  desired  shape  of  the  free  surface 

in  the  final  product  (Example  5). 


Figure  24.  The  objective  function  and  the  norm  of  the  gradient  of  the  objective  function  versus  the  number 
of  iterations  for  the  open-forging  die  design  problem  (Example  5). 
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5.  CONCLUSIONS 

An  accurate  and  efficient  continuum  sensitivity  analysis  was  presented  for  the  design  of  metal 
forming  processes.  The  developed  methodology  is  applicable  for  the  calculation  of  the  sensitivities 
of  the  deformation,  stresses  and  state  variables  with  respect  to  variations  in  the  die  surface  as  well 
as  variations  in  other  (non-shape  related)  process  parameters. 

An  effective  method  was  presented  for  computing  the  sensitivities  of  the  contact  tractions.  The 
method  is  developed  by  differentiation  of  the  continuum  contact  and  frictional  constraints  after  a 
number  of  regularizing  assumptions  were  introduced.  The  linearity  of  the  sensitivity  algorithm  was 
preserved  and  no  iterations  or  augmentations  were  required  in  the  calculation  of  various  sensitivity 
fields. 

Various  examples  were  presented  to  demonstrate  the  accuracy  and  potential  of  the  developed 
sensitivity  analysis.  Both  extrusion  and  open-die  forging  processes  were  considered.  Preform  design 
problems  based  on  the  present  analysis  are  addressed  in  a  companion  paper  [33], 


APPENDIX  A.  DERIVATION  OF  THE  LINEAR  RELATION  BETWEEN  ft+1  AND  F„+i 

*  * 

A  sequence  of  steps  is  followed  in  order  to  express  F^+1  as  a  linear  function  of  F„+i.  The  variable 
D„+i  is  introduced  as  follows: 


d«+i  =  [n;+,(F«+.  f„+')f,:+i 


(Al) 


and  G„+i  is  defined  as  follows: 

G„+i  =  D„+]  —  C„+i  (A2) 

where  C„+i  is  given  by  Equation  (36).  For  an  isotropic  hyperelastic  model  with  shear  modulus 
"V,  one  can  write  the  following: 

*  *  * 

T;,+i  =  2<men+l  -  3  tr(Ee„+l  )I  (A3 ) 

As  shown  in  Reference  [26], 

tr(fe+1)  =  tr(U„e+i[Ue]„-+11)  (A4) 

Using  the  tensor  property  tr(AB)  =  tr(BA)  and  the  sensitivity  of  the  the  polar  decomposition  of 
F)h  I,  one  can  show  that 


tr(F„e+l[Fe];+i )  =  tr(U„e+1[Ue]-11 ) 
and  using  Equation  (A4)  conclude  that 

tr(F„e+1[Fe]„-+11)  =  tr(le„+1) 


(A5) 


(A6) 
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Using  Equations  (41),  (A2)  and  (A6)  and  the  fact  that  deviatoric  tensors  are  traceless,  one  arrives 
at  the  following  equation: 


tr  (E® +i )  =  tr(G„+i )  (A7 ) 

* 

By  considering  the  dot  product  of  Equation  (41)  with  T'+1,  one  can  eliminate  <?„+ 1  from  this 
equation.  This  is  shown  by  the  sequence  of  steps  that  follow. 

-i  * 

Let  us  first  consider  the  quantity,  ([Fe]H+]F^+i )  •  T'i+1.  Using  the  sensitivity  of  the  polar  decom¬ 
position  of  F®+1,  one  can  write, 


{[Fe]^,Fne+,}  •  f;,+l  =  {[ux+,([Re]I+,R«+1)u„e+1}  •  rn+]  +  ([Ue];+,u®+,)-  rn+]  (A8) 


Since  U®+1  and  E®+1  commute  in  multiplication,  U®+1  and  T„+i  and  hence  U®+1  and  T;'+ ,  also 
commute  in  multiplication.  Thus, 


{[Ue]J^(11([Re]J_l_IR%+i  )U„e+1 }  •  T'n+l 

=  {[Re]„T+1Re„+i}  •  {[ue]^,f:+lu;+1} 

=  {[Re]„T+1Re„+1}-t;1+1 

=  0  (A9) 

where  the  skew  symmetry  of  [Re]„+1R^+i  and  symmetry  of  T'+1  were  taken  into  account.  The 
following  simplified  equation  is  finally  obtained  from  Equation  (A8), 

([Fe];+,F«+| )  •  T;i+1  =  ([Ue]^,U®+1 )  .  T'n+l  (A10) 

Using  the  diagonal  representation  of  the  symmetric  tensor  U®+1,  one  can  write 

U;+,=Q„+|A„+|QI+,  (All) 

where  A„+i  is  a  diagonal  tensor  and  Q„+i  is  a  rotation  tensor.  Therefore, 

U®+1  =  Q„+iA„+iQj+1  +  2  sym(Q„+1  Ql+iU„e+1 )  (A12) 

Since  U,®+1  and  T'+]  commute  in  multiplication  and  both  are  symmetric,  (U®+1)_1T';+1  is 
* 

symmetric.  Since  Q„+iQ„+i  is  skew  symmetric  and  T';+1  symmetric,  the  following  equation  holds: 

{[UeC  sym(Q„+1  Qj+1U„e+1)}  •  %+1  =  0  (A13) 

Thus  using  Equations  (A12)  and  (A13),  one  can  write 

([Ue]^,U*+1 ) .  t;(+l  =  (Q„+1  A„+1  A-\ Ql+I  )-  f'+1  (A14) 
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Similarly,  one  can  write  the  following: 


E^+i  =  ln(U*+i )  =  Qn+i  ln(A„+1  )Q 


«+ 1 


E«+i  —  Q«+iA„+iA„+1Q„+1  +  2  sym(Q„+1  Qn+1En+1) 


(A15) 

(A16) 


with 


sym(Q„+1  QI+1E:+i)  •  f;i+l  =  0 


Therefore  using  Equations  (A16)  and  (A17),  one  derives  the  following: 


(A17) 


E  n+i  •  Tn+1  —  (Q„+iA„+iA„+1Qn+1)  •  T„+] 
Combining  Equations  (A14)  and  (A18),  one  can  write  the  following: 

*  * 

Also  following  Equation  (25)  and  using  T'„+i  =2@Ee'n+1,  one  can  write 


* 
<T/i+ 1 


* 

Ec/ 


M+l’ 

Y4  * 

T'  _  J  J  l?e  T' 

n+ 1"  1  n+l  —  ~  J1'  n+i"  1h+1 


<7„+i  n„+i 

The  following  simplification  is  obtained  using  Equations  (A10),  (A19)  and  (A20): 

3^ 


([Fe],7+|E«+|)‘  t;1+1  =  k+ r  t;+1  =  ^<T„+1 


(A18) 


(A19) 


(A20) 


(A21) 


Taking  the  dot  product  of  Equation  (A22)  with  T')+1  and  using  Equations  (14),  (25)  and  (A21), 
results  in  the  following: 


~  J  <7«+i  2  <T„+ifl„+i  2dH_|_jb„+i 

(j«+r  l„+i  =  o-»+i  <  H - 3 - H - 3 - 


(A22) 


and  thus  an+l  is  finally  expressed  as 


* 

®n+ 1 


3^G„+1-T'+1 


<b,+i (1  T  2i^(nw+i  T  <?„+ib„+i)) 


(A23) 


Let  us  now  define  the  following  tensor  (function  of  F,®+i ): 

EI„+i  =  [Fe]„+11F^+i  +  2^a„+iE‘;,+i 


(A24) 


Using  this  definition  and  Equations  (A2),  (A3),  (A7),  Equation  (41)  can  now  be  simplified  as 
follows: 


2^a„+itr(G„+i) 

*1„+1  —  vr„+i  -\ - - - 1  ■ 


Y$bn+\G, 


ln+l'J«+l  ‘  F)i+] 


<7„+i(l  T  2^(a„+i  T<7«+ibw+i)) 


T 


«+ 1 


(A25) 
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Table  Al.  Material  parameters  for  A1  1100-0 
at  an  initial  temperature  673  K. 


Material  parameter 

Value 

A 

1.90593  x  107  s-1 

Q/R 

2.1086  x  104K“‘ 

7.0 

m 

0.23348 

so 

29.7  MPa 

ho 

1115.6  MPa 

a 

1.3 

s 

18.92  MPa 

n 

0.07049 

20.2  GPa 

Jf 

66.0  GPa 

or  finally  as 


H„+i  —  Dw-|_i  -\~ 


2 ^an+l  tr(D„+1) 


I 


YSbn+\  D, 


n+l*-*n+l 


T' 

H+l 


@n+  l(l  T"  2 A( il n  |  l  -\~  6n-Y\bn+\  )) 


T'„ 


n+ 1 


—  C„+1  + 


3%+iQ+rT; 


/z+l 


^n-t-iO  T-  2Y/(c/„  |  ]  -(-  )) 


T' 


n+1 


(A26) 


Notice  that  D„+i  is  linearly  related  with  Fn+i  and  C„+i  is  independent  of  Fn+i-  With  the  definition 
of  H„+i  from  Equation  (A24),  one  can  transform  Equation  (A26)  to  the  form  of  Equation  (42). 

Also,  given  F„+i  and  using  the  definition  of  D„+i  and  C„+i  from  Equations  (Al)  and  (36), 

* 

Equation  (42)  can  be  solved  numerically  to  obtain  F,®+i. 


APPENDIX  B.  CONSTITUTIVE  MODEL  FOR  Al-1100  AT  673  K 


The  constitutive  model  for  1100-A1  at  0  —  673  K  is  here  taken  from  [27],  The  flow  function  /  is 
given  as  follows: 


/(( j,  s,  6)  =  A  exp 


and  the  hardening  function  g  is  given  by 

g{a ,  s,  9)  =  h(o,s)f  {a,  s,  9) 
where  the  function  h  is  defined  as  follows: 


(Bl) 


(B2) 


/z(<7,s)  =  /*o|l  — 


(B3) 
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with 


s  =s 


f(o,s,  6) 


exp 


n 


The  specific  values  of  the  mechanical  and  thermal  parameters  are  given  in  Table  Al. 


(B4) 
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